'''
Created on Oct 16, 2009

@author: oabalbin
'''
import operator
from operator import itemgetter
import sys
import glob
import os
import math
import copy
import pickle
import heapq
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import numpy as np
import scipy as sp
from scipy import stats
from collections import deque, defaultdict

from common.classes import myArray, gene_copa
import signatures.db.gene_annotations as ga


class ExpressionData(object):
    '''
    classdocs
    '''


class Sample:
    '''
    
    '''
    def __init__(self,  mysamplename, mygene_raw_cov, mygenelist, mygene_normalized_cov, my_sample_norm_covDic):
        #self.gene_rawcov = mygene_raw_cov
        self.geneRank = []
        self.mygenes = mygenelist
        self.samplename = mysamplename
        # raw values
        self.gene_rawcov = mygene_raw_cov
        # normalization accross samples
        self.gene_norcov = mygene_normalized_cov
        # normalization in the sample
        self.gene_sample_normcov = my_sample_norm_covDic
        # Expression of all genes in sample
        self.mygeneDic = self.make_mygeneDic()
        self.myrawgeneDic = self.make_my_rawgeneDic()
        self.mysample_genecovDic = self.make_my_samplegeneDic()
        self.outliers= {}
                
        self.vector_of_outliers = np.empty(len(mygenelist))
        
    def make_mygeneDic(self):
        mygeneDic={}
        for gene, expr in zip(self.mygenes,self.gene_norcov[self.samplename]):
            #if self.samplename == "HCC1954":
                #print gene, expr 
            mygeneDic[gene] = expr 
        return mygeneDic
    
    def  make_my_rawgeneDic(self):
        mygeneDic={}
        for gene, expr in zip(self.mygenes,self.gene_rawcov[self.samplename]):
            mygeneDic[gene] = expr     
        return mygeneDic
    
    def  make_my_samplegeneDic(self):
        mygeneDic={}
        for gene, expr in zip(self.mygenes,self.gene_sample_normcov[self.samplename]):
            mygeneDic[gene] = expr     
        return mygeneDic
    

class GeneArray:
    '''
    
    '''
    def __init__(self,  expvalues, geneValueInSample, geneZscoreSample, genemedian, genesd, celllines, genelist, sampleCrudeExpValues=[]):
        
        self.geneRawValues = expvalues
        self.geneExpInSample = geneValueInSample
        self.geneZscoreInSample = geneZscoreSample
        self.genemed = genemedian
        self.genesdt = genesd
        self.celllines = celllines
        self.genelist = genelist
        self.sampleCrudeExpValues=sampleCrudeExpValues
    


def open_files(folderpath,filename):
    '''
    Open the file with the expression for the kinases and other genes.
    '''
    myfiles=[]
    for infile in glob.glob( os.path.join(folderpath, filename) ):
        myfiles.append(infile)
    
    return myfiles


def create_mysamples(myfiles,normal_samples=[]):
    '''
    It creates a sample object and calculates a z-score for the expression level of each gene
    in the sample. Z-score (expvalues - genemedian)/genesd
    Returns a geneArray and also a samplelist with the expression for each sample
    The input matrix is already lo2 normalized
    '''
    geneRawValues,gene_norm_cov, geneZscoreInSample, genemed, genesdt, celllines = {},{},{},{},{}, deque()
    
    
    for thisfile in myfiles:
        genelist=[]
        file = open(thisfile,'r')
        first=True

        for line  in file:

            if line is None:
                continue 
            if line.startswith('#'):
                continue
            if line.startswith('gene'):
                celllines = line.strip('\n').split('\t')[2:]
                if normal_samples:
                    nsp_index=[]
                    for tsp in normal_samples:
                        print tsp
                        nsp_index.append(celllines.index(tsp))
                    #print tsp
                    #[nsp_index.append(celllines.index(tsp)) for tsp in normal_samples]
                    #    nsp_index.append(celllines.index(sp))
                continue
            
            fields = line.strip('\n').split('\t') 
            
            genename = fields[0]
            genelist.append(genename)
            
            exp=np.zeros(len(fields[2:]))
            expmasked = np.ma.make_mask_none((len(fields[2:]),))
            for i, e in enumerate(fields[2:]):
                if e=="--" :
                    expmasked[i] = False
                    exp[i] = float(0.1)
                else:
                    exp[i] = float(e) + 0.1
                    #exp.append(float(i))
            #exp = map(float,fields[2:])
            # shift the values by one. In order to avoid the zero coverage values which will be inf when log is applied.
            #expvalues = np.log2(np.array([exp]) + 2)
            expvalues = np.ma.array([exp],mask=[expmasked])
            expvalues = np.array([exp])
            # Calculating gen median and std
            normals_mean = np.ma.mean(expvalues)/np.ma.mean(expvalues[:,nsp_index])
            #expvalues = expvalues/normals_mean
             
            genemean = np.ma.mean(expvalues)
            genemedian = np.ma.median(expvalues)
            genesd = np.ma.std(expvalues)
            uppercent=75
            gene75percent = sp.stats.scoreatpercentile(expvalues[0,:],uppercent)

            # Calculating scores
            # you are using log2 values
            geneValueInSample = expvalues - genemedian
            #expvalues = expvalues - genemedian
            expvalues = expvalues - gene75percent
            #geneZscoreSample= expvalues - genemedian
            geneZscoreSample = (expvalues - genemedian)/genesd
            #geneZscoreSample = (expvalues - genemean)/genesd
            #print np.mean(geneZscoreSample), np.std(geneZscoreSample)
            
            #Creating class geneArray
            geneRawValues[genename] = expvalues                      # normalized over the expression value of the normals
            gene_norm_cov[genename] = geneValueInSample
            geneZscoreInSample[genename] = geneZscoreSample
                       
            genemed[genename] = genemedian
            genesdt[genename] = genesd
            
            
            #Creating an array of expression organized by sample
            # This could be better organized into a numpy array of defined dimensions.
            if first==True:
                sampleRawValues = np.array(expvalues)
                #samplegene_norm_cov = np.array(geneValueInSample)
                samplegene_norm_cov = np.array(geneZscoreSample)
                first=False
            else:
                sampleRawValues = np.concatenate((sampleRawValues,expvalues),axis=0)
                #samplegene_norm_cov = np.concatenate((samplegene_norm_cov,geneValueInSample),axis=0)
                samplegene_norm_cov = np.concatenate((samplegene_norm_cov,geneZscoreSample),axis=0)
    
    # Create an GeneArray object
    thisGeneArray = GeneArray(geneRawValues, gene_norm_cov, geneZscoreInSample, genemed, genesdt, celllines, genelist)
    # Create a list of celllines
    
    
    samplelist=[]
    mysamples_array= np.transpose(np.zeros((len(genelist) ,len(celllines))))
     
    for i,cell in enumerate(celllines):
        cellind = celllines.index(cell)
        geneDic, gene_norm_covDic, sample_norm_covDic= defaultdict(),defaultdict(),defaultdict()
        
        thissample_raw_values = sampleRawValues[:,cellind]
        #thissample_raw_values = samplegene_norm_cov[:,cellind]
        geneDic[cell] = thissample_raw_values #sampleRawValues[:,cellind]
        gene_norm_covDic[cell] = samplegene_norm_cov[:,cellind]
        
        #normalization of gene expression within a particular sample. 
        gene_median = np.median(thissample_raw_values)
        gene_mad = np.median(np.abs(thissample_raw_values - gene_median))   
        gene_std = np.std(thissample_raw_values)
        gene_insample_normalization = (thissample_raw_values - gene_median)/gene_mad
        #gene_insample_normalization = (thissample_raw_values - gene_median)/gene_mad
        
                
        sample_norm_covDic[cell] = gene_insample_normalization
        mysamples_array[i,:] = gene_insample_normalization
        
        thissample = Sample(cell, geneDic, genelist, gene_norm_covDic, sample_norm_covDic)
        samplelist.append(thissample)
    
    return thisGeneArray, samplelist, np.transpose(mysamples_array)
            

def outliers(samplelist, uppercent, celllines, genelist, norm_accross_sample=False):
    '''
    It creates a list of candidates based on the zscoreof each gene.
    The candidate list is constituted by the genes with zscore >= uppercent
    '''
    
    # matrix for reporting the presence or absence of a particular outlier
    outlier_matrix = np.zeros((len(genelist),len(celllines)))
    
    for thissample in samplelist:
        '''
        if thissample.samplename=="HCC1428":
            exprlist = thissample.gene_sample_normcov[thissample.samplename]
            print np.sort(exprlist)
        ''' 
        
        if norm_accross_sample:
            # It uses the normalization accros samples but for each gene. zscore for each genes.
            exprlist = thissample.gene_norcov[thissample.samplename]
            my_expression_dictionary = thissample.mygeneDic
        
        else:
            # It uses the normalization inside the samples
            exprlist = thissample.gene_sample_normcov[thissample.samplename]
            my_expression_dictionary = thissample.mysample_genecovDic
        
        #uppercentscore = sp.stats.scoreatpercentile(abs(exprlist),uppercent)
        uppercentscore = sp.stats.scoreatpercentile(exprlist,uppercent)
        #downperexonscore = sp.stats.scoreatpercentile(a_exprlist,downpercent)
        #uppercentscore2 = sp.stats.scoreatpercentile(thissample.mygeneDic.values(),uppercent) 
        
        candidates = {}
        #for gene, expr_val in thissample.mygeneDic.iteritems():
        for gene, expr_val in my_expression_dictionary.iteritems():      
            #print gene, expr_val, uppercentscore
            #if (abs(expr_val) >= uppercentscore):
            if (expr_val >= uppercentscore):
                #print "accepted value "
                #print expr_val
                candidates[gene] = expr_val
            else:
                continue
    
        #sort candidates
        #print sorted(candidates.items(), key=operator.itemgetter(1), reverse=True)
        print thissample.samplename, sorted(candidates.items(), key=operator.itemgetter(1), reverse=True)
        
        # Assign the candidates to the gen outliers in the sample.
        thissample.outliers = candidates
        
        ### Creating a matrix of outliers-sample of presence or abscence of that particular gene
        
        for gen, expr_val in candidates.iteritems():
            i = genelist.index(gen)
            j = celllines.index(thissample.samplename)
            outlier_matrix[i,j] = expr_val #thissample.myrawgeneDic[gene] #expr_val
        
    return outlier_matrix
    

def sample_group_by_outlier_gene(samplelist, gene):
    '''
    It generates a group of samples based on the presence of a particular gene
    in the set of candidate genes of each sample.
    returns a list with the group of samples.
    '''
    
    sample_group = [thissample.samplename for thissample in samplelist if gene in thissample.outliers.keys()]
    
    return sample_group
    

def recurrent_set_of_genes(samplelist):
    '''
    It find the most outlier expressed genes in each sample. Then it looks for the most common outlier gene accross 
    samples. Iterate the process until there is not recurent outlier gene. 
    '''
    
    '''
    Zero:
    for each sample get the list of candidates:
    first_ouliers = list(get pair  (gene,value) the maximum outlier of each sample)
    determine the more recurrent gene among first_ouliers
    Discard from next analysis samples in which that gene is not the most outlier, 
    Discard totally sample in which that gene is not in the candidates sample
    
    Option One:
    For samples with outlier:
    Intersect the list of remainder candidates genes to find a common set of "affected" genes
     
    Option Two:
    Repeat Zero with the second most outlier gene.
    until not recurrent outliers can be found
    '''
    outliers_list={}
    sorted_candidates={}
    samplelist_not_outliers=[]
    for thissample in samplelist:
        sorted_candidates[thissample] = sorted(thissample.outliers.items(), key=operator.itemgetter(1), reverse=True)
        
        if len(sorted_candidates[thissample]) > 0:
            gene = (sorted_candidates[thissample])[0][0]
            
            if gene in outliers_list:
                recurrency = outliers_list.get(gene)
                outliers_list[gene] = recurrency+1
            else:
                outliers_list[gene] = 1
        else:
            samplelist_not_outliers.append(thissample.samplename)
             
    #print outliers_list
        
    recurrent_list = heapq.nlargest(1, outliers_list.iteritems(), operator.itemgetter(1))
    
    print  recurrent_list

    
def copa(gene_covs):
    '''
    Calculates the copa index in order to make an outlier analysis. 
    '''
    gene_median = np.median(gene_covs,axis=0)    
    gene_mad = np.median(np.abs(gene_covs - gene_median))
    gene_copa = (gene_covs - gene_median) / gene_mad
    
    return gene_copa, gene_median



def plot_outliers(samplelist,geneArray, folderpath):
    '''
    for top 10 expressed Kinases in a particular samples plot 
    expression value for that gene and the median expression value
    across the samples
    '''
    sorted_candidates={}
    
    for thissample in samplelist:
        sorted_candidates[thissample.samplename] = deque(sorted(thissample.outliers.items(), key=operator.itemgetter(1), reverse=True))
    
        #exp_in_sample, exp_in_set, std_in_set = {},{},{}
        noutliers=len(sorted_candidates[thissample.samplename])       
        exp_in_sample, exp_in_set, std_in_set, thisgenes = np.empty(noutliers),np.empty(noutliers),np.empty(noutliers),[]
        
        if len(sorted_candidates[thissample.samplename]) > 0:
            for i in range(noutliers):# sorted_candidates[thissample.samplename]:
                #print genepair
                genepair = sorted_candidates[thissample.samplename][i]
                gene = (genepair)[0]
                thisgenes.insert(i,gene)
                exp_in_sample[i] = (genepair)[1]
                exp_in_set[i] = geneArray.genemed[gene]
                std_in_set[i] = geneArray.genesdt[gene]
                '''
                gene = (genepair)[0]
                # To plot the z-score or any normalized value
                exp_in_sample[gene] = (genepair)[1] 
                # To plot the expression value of the gen in the sample. 
                #exp_in_sample[gene] = thissample.myrawgeneDic[gene]
                exp_in_set[gene] = thissample.myrawgeneDic[gene]
                #exp_in_set[gene] = geneArray.genemed[gene]
                #std_in_set[gene] = geneArray.genesdt[gene]
                std_in_set[gene] = 0.0
                #print exp_in_set[gene],  genemedian
                '''
                
            plot_sample_centric_view(exp_in_sample, exp_in_set, std_in_set, thissample.samplename, thisgenes, folderpath) 
    

def plot_sample_centric_view(exp_in_sample, exp_in_set, std_in_set, sample, genenames, folderpath):
    
 
    ind = np.arange(len(exp_in_sample))  # the x locations for the groups
    width = 0.45                         # the width of the bars
    #print std_in_set.values() #, ind,exp_in_sample.values()
    #print len(exp_in_set.values()), len(std_in_set.values())

    plt.bar(ind,exp_in_sample,width,color='r')  
    #plt.aspect(1.)
    #plt.bar(ind+width,exp_in_set,width, color='b',yerr =std_in_set)#, yerr =std_in_set.values())
    #plt.xlabel('Genes')
    plt.ylabel('Z-score (red), Log2 Expression (blue)')
    #axScatter.set_xlim(xlim)
    #axScatter.set_ylim(ylim) 
    plt.title(sample)
    plt.xticks(ind+width, genenames, rotation='vertical', size=6)
    #axScatter.set_xticklabels( exp_in_set.keys(), rotation='vertical', size=6)

    '''
    ax1 = plt.subplot(111)
    ax1.bar(ind,exp_in_sample,width,color='r')  
    #plt.aspect(1.)
    ax1.set_ylabel('Z score')
    ax2 = plt.twinx()
    ax2.bar(ind+width,exp_in_set,width, color='b',yerr =std_in_set)#, yerr =std_in_set.values())
    #plt.xlabel('Genes')
    ax2.set_ylabel('Log2 Expression')
    ax2.yaxis.tick_right()
    #axScatter.set_xlim(xlim)
    #axScatter.set_ylim(ylim) 
    ax2.set_title(sample)
    ax2.set_xticks(ind+width, genenames)#, rotation='vertical', size=6)
    #axScatter.set_xticklabels( exp_in_set.keys(), rotation='vertical', size=6)
    plt.setp( ax2.get_xticklabels(), rotation = 'vertical', fontsize =6 )

    
    axScatter = plt.subplot(111)
    axScatter.bar(ind,exp_in_sample.values(),width,color='r')  
    axScatter.set_aspect(1.)
    axScatter.bar(ind+width,exp_in_set.values(),width, color='b',yerr =std_in_set.values())#, yerr =std_in_set.values())
    axScatter.set_xlabel('Genes')
    axScatter.set_ylabel('Log2 Expression')
    #axScatter.set_xlim(xlim)
    #axScatter.set_ylim(ylim) 
    axScatter.set_title(sample)
    axScatter.set_xticks(ind+width)
    axScatter.set_xticklabels( exp_in_set.keys(), rotation='vertical', size=6)
    '''
    
    plt.savefig(folderpath+sample+'_3.png')
    plt.close()


def write_dict2file(outfile,thisDict):
    """
    """
    for mykey, myvals  in thisDict.iteritems():
        mykey = str(mykey)
        if myvals.shape[1] != 1:
            el = map(str,myvals[0,:])
            outfile.write(mykey+'\t'+",".join(el).replace(',','\t')+'\n')
        else:
            el = str(myvals)
            outfile.write(mykey+'\t'+str(myvals)+'\t'+'\n')


def write_dict2file2(outfile,thisDict,select_names):
    """
    """
    for mykey, myvals  in thisDict.iteritems():
        mykey = str(mykey)
        
        if mykey in select_names.keys():
            load = select_names[mykey]
            if myvals.shape[1] != 1:
                el = map(str,myvals[0,:])
                outfile.write(mykey+'\t'+load+'\t'+",".join(el).replace(',','\t')+'\n')
            else:
                el = str(myvals)
                outfile.write(mykey+'\t'+str(myvals)+'\t'+'\n')


def write_array2file(outfile, H, rownames=[], select_names=[]):
    """
    It writes the numpy array to a file. 
    """
    
    for i in range(H.shape[0]):
        el = map(str,list(H[i,:]))
        if rownames:
            if select_names:
                if rownames[i] in select_names:
                    outfile.write('T'+'\t'+str(rownames[i])+"\t"+",".join(el).replace(',','\t')+'\n')
                else:
                    outfile.write('N'+'\t'+str(rownames[i])+"\t"+",".join(el).replace(',','\t')+'\n')
            else:
                #print str(rownames[i])+"\t"+",".join(el).replace(',','\t')
                outfile.write(str(rownames[i])+"\t"+",".join(el).replace(',','\t')+'\n')
        else:
            outfile.write(",".join(el).replace(',','\t')+'\n')
            

def write_array2file2(outfile, H, rownames=[], select_names={}):
    """
    It writes the numpy array to a file. 
    """
    
    for i in range(H.shape[0]):
        el = map(str,list(H[i,:]))
        if rownames:
            if select_names:
                if rownames[i] in select_names.keys():
                    #print str(rownames[i])+"\t"+",".join(el).replace(',','\t')
                    load = select_names[rownames[i]]
                    outfile.write(str(rownames[i])+"\t"+load+"\t"+",".join(el).replace(',','\t')+'\n')
        else:
            outfile.write(",".join(el).replace(',','\t')+'\n')


def list_of_names_dict(inputfile):
        
        genlist=defaultdict()
        for line in inputfile:        
            line = line.strip('\n')
            fields = line.split('\t')
            #print fields 
            # To skip headers, star reading samples in column 7 of the file
            if fields[0][0] == '#': 
                continue                             
            if fields[0]=='':
                continue
            else:
                genlist[fields[0]]=fields[1]
        
        return genlist

def list_of_names(inputfile):
        
        genlist=deque()
        for line in inputfile:        
            line = line.strip('\n')
            fields = line.split('\t')
            #print fields 
            # To skip headers, star reading samples in column 7 of the file
            if fields[0][0] == '#': 
                continue                             
            if fields[0]=='':
                continue
            
            genlist.append(fields[0])
        
        return genlist



def print_sample_normalized_matrix(list_of_samples,cellsnames,geneslist,outfile,select_names={}):
    """
    """
    
    mysamples_array= np.zeros((len(list_of_samples),len(geneslist)))
    sample_expression_dictionary = defaultdict(list)
    for i,thissample in enumerate(list_of_samples):
        for gen in geneslist:
            sample_expression_dictionary[gen].append(thissample.mysample_genecovDic[gen])
                
    if select_names:
        outfile.write("geneName"+'\t'+"load"+'\t'+",".join(cellsnames).replace(',','\t')+'\n')
    else:
        outfile.write("geneName"+'\t'+",".join(cellsnames).replace(',','\t')+'\n')
    
    for mykey, myvals  in sample_expression_dictionary.iteritems():
        if select_names:
            if mykey in select_names.keys():
                load = select_names[mykey]
                el = map(str,myvals)
                outfile.write(mykey+'\t'+load+'\t'+",".join(el).replace(',','\t')+'\n')
        else:
            el = map(str,myvals)
            outfile.write(mykey+'\t'+",".join(el).replace(',','\t')+'\n')




def outliers_atpercentile(samplelist, uppercent, celllines, genelist, outlier_matrix, norm_accross_sample=False):
    '''
    It creates a list of candidates based on the zscoreof each gene.
    The candidate list is constituted by the genes with zscore >= uppercent
    '''
    
    # matrix for reporting the presence or absence of a particular outlier
    #outlier_matrix = np.zeros((len(genelist),len(celllines)))
    
    for thissample in samplelist:
        
        if norm_accross_sample:
            # It uses the normalization accros samples but for each gene. zscore for each genes.
            exprlist = thissample.gene_norcov[thissample.samplename]
            my_expression_dictionary = thissample.mygeneDic
        
        else:
            # It uses the normalization inside the samples
            exprlist = thissample.gene_sample_normcov[thissample.samplename]
            my_expression_dictionary = thissample.mysample_genecovDic
        
        #uppercentscore = sp.stats.scoreatpercentile(abs(exprlist),uppercent)
        uppercentscore = sp.stats.scoreatpercentile(exprlist,uppercent)
        #downperexonscore = sp.stats.scoreatpercentile(a_exprlist,downpercent)
        #uppercentscore2 = sp.stats.scoreatpercentile(thissample.mygeneDic.values(),uppercent) 
        
        candidates = {}
        #for gene, expr_val in thissample.mygeneDic.iteritems():
        for gene, expr_val in my_expression_dictionary.iteritems():      
            #print gene, expr_val, uppercentscore
            #if (abs(expr_val) >= uppercentscore):
            if (expr_val >= uppercentscore):
                #print "accepted value "
                #print expr_val
                candidates[gene] = expr_val
            else:
                continue
    
        sorted_candidates = sorted(candidates.items(), key=operator.itemgetter(1), reverse=True)
        # Assign the candidates to the gen outliers in the sample.
        thissample.outliers = candidates
        
        ### Creating a matrix of outliers-sample of presence or abscence of that particular gene
        outlier_genes=[]
        j = celllines.index(thissample.samplename)
        for rank, pair in enumerate(sorted_candidates):
            i = genelist.index(pair[0])
            expval=pair[1]
            outlier_genes.append(i)
            #j = celllines.index(thissample.samplename)
            outlier_matrix[i,j] += expval
        
        outlier_genes = np.array(outlier_genes)
        #outlier_matrix[-outlier_genes,j] += 100 
        '''
        if thissample.samplename == "ZR-75-30":
            genelist2 = np.array(genelist)
            print genelist2[outlier_genes]
            print outlier_matrix[outlier_genes,j]
        '''
        
        
    return outlier_matrix





def create_mysamples_atpercentile(myfiles, percentile, normal_samples=[]):
    '''
    It creates a sample object and calculates a z-score for the expression level of each gene
    in the sample. Z-score (expvalues - genemedian)/genesd
    Returns a geneArray and also a samplelist with the expression for each sample
    The input matrix is already lo2 normalized
    
    Note: jun-4-2010. 
    This functions sguld be divided into 2 different ones. The first should generate the GeneArray structure and another structure
    with the sampleRawValues,  samplegene_norm_cov, sampleCrudeExpValues. Similar to MyArray structures so I can operate on that array later.
    The second function shoul mahe the sample normalization, and the copa outlier inside the sample and accross the sample. 
    '''
    geneRawValues,gene_norm_cov, geneZscoreInSample, crudeRawValues, genemed, genesdt, celllines = {},{},{},{},{},{}, deque()
    
    
    for thisfile in myfiles:
        genelist=[]
        file = open(thisfile,'r')
        first=True

        for line  in file:

            if line is None:
                continue 
            if line.startswith('#'):
                continue
            if line.startswith('gene'):
                celllines = line.strip('\n').split('\t')[2:]
                if normal_samples:
                    nsp_index=[]
                    [nsp_index.append(celllines.index(tsp)) for tsp in normal_samples]
                    #    nsp_index.append(celllines.index(sp))
                continue
            
            fields = line.strip('\n').split('\t') 
            
            genename = fields[0]
            genelist.append(genename)
            
            exp=np.zeros(len(fields[2:]))
            expmasked = np.ma.make_mask_none((len(fields[2:]),))
            for i, e in enumerate(fields[2:]):
                if e=="--" :
                    expmasked[i] = False
                    #exp[i] = float(0.1)
                    exp[i] = float(2.0)
                    #exp[i] = np.log2(float(0.1))
                else:
                    #exp[i] = float(e) + 0.1
                    exp[i] = float(e) + 2.0
                    #exp.append(float(i))
            #exp = map(float,fields[2:])
            # shift the values by one. In order to avoid the zero coverage values which will be inf when log is applied.
            #expvalues = np.array([exp])
            expvalues = np.log2(np.array([exp]))
            expvalues = np.ma.array(expvalues,mask=[expmasked])
            
            #expvalues_crude=expvalues
            
            # Calculating gen median and std
            normals_mean = np.ma.mean(expvalues)/np.ma.mean(expvalues[:,nsp_index])
            #expvalues = expvalues/normals_mean 
            genemean = np.ma.mean(expvalues)
            genemedian = np.ma.median(expvalues)
            genesd = np.ma.std(expvalues)

            gene_percentile = sp.stats.scoreatpercentile(expvalues[0,:],percentile)

            # Calculating scores
            # you are using log2 values
            # For log2 setup
            expvalues_crude= expvalues #- genemedian
            gene_mad_pop = np.median(np.abs(expvalues - genemedian))
            expvalues_perc = expvalues - gene_percentile #)/gene_mad_pop
            # For raw values
            #expvalues_crude=expvalues#/genemedian
            #geneValueInSample = expvalues/genemedian
            #expvalues_perc = expvalues/gene_percentile
            #geneZscoreSample= expvalues - genemedian
            geneValueInSample = expvalues - genemedian
            geneZscoreSample = (expvalues - genemedian)/genesd
            #geneZscoreSample = (expvalues - genemean)/genesd
            #print np.mean(geneZscoreSample), np.std(geneZscoreSample)
            
            # Raw raw data, not normalized
            crudeRawValues[genename]=expvalues_crude
            #Creating class geneArray
            geneRawValues[genename] = expvalues_perc                      # normalized over the expression value of the normals
            gene_norm_cov[genename] = geneValueInSample
            geneZscoreInSample[genename] = geneZscoreSample
                       
            genemed[genename] = genemedian
            genesdt[genename] = genesd
            
            
            #Creating an array of expression organized by sample
            # This could be better organized into a numpy array of defined dimensions.
            if first==True:
                sampleRawValues = np.array(expvalues_perc)
                #samplegene_norm_cov = np.array(geneValueInSample)
                samplegene_norm_cov = np.array(geneZscoreSample)
                sampleCrudeExpValues = np.array(expvalues_crude)
                first=False
            else:
                sampleRawValues = np.concatenate((sampleRawValues,expvalues_perc),axis=0)
                #samplegene_norm_cov = np.concatenate((samplegene_norm_cov,geneValueInSample),axis=0)
                samplegene_norm_cov = np.concatenate((samplegene_norm_cov,geneZscoreSample),axis=0)
                # Crude expression in sample
                sampleCrudeExpValues = np.concatenate((sampleCrudeExpValues,expvalues_crude),axis=0)
    
    # Create an GeneArray object
    thisGeneArray = GeneArray(geneRawValues, gene_norm_cov, geneZscoreInSample, genemed, genesdt, celllines, genelist, sampleCrudeExpValues)
    # Create a list of celllines
    
    
    ######## In sample processing.
        
    ### Include for the analysis of the outlier only genes whose median is above the 25 percentile of the medians of all other genes.
    med_percentile=25
    median_percentile = sp.stats.scoreatpercentile(genemed.values(),med_percentile)
    
    samplelist=[]
    mysamples_array= np.transpose(np.zeros((len(genelist) ,len(celllines))))
     
    for i,cell in enumerate(celllines):
        cellind = celllines.index(cell)
        geneDic, gene_norm_covDic, sample_norm_covDic= defaultdict(),defaultdict(),defaultdict()
        
        thissample_raw_values = sampleRawValues[:,cellind]
        #thissample_raw_values = samplegene_norm_cov[:,cellind]
        geneDic[cell] = thissample_raw_values #sampleRawValues[:,cellind]
        gene_norm_covDic[cell] = samplegene_norm_cov[:,cellind]
        
        #normalization of gene expression within a particular sample. 
        gene_median = np.median(thissample_raw_values)
        gene_mad = np.median(np.abs(thissample_raw_values - gene_median))   
        gene_std = np.std(thissample_raw_values)
        gene_insample_normalization = (thissample_raw_values - gene_median)/gene_mad
        #gene_insample_normalization = (thissample_raw_values - gene_median)/gene_mad
        
                
        sample_norm_covDic[cell] = gene_insample_normalization
        mysamples_array[i,:] = gene_insample_normalization
        
        thissample = Sample(cell, geneDic, genelist, gene_norm_covDic, sample_norm_covDic)
        samplelist.append(thissample)
    
    return thisGeneArray, samplelist, np.transpose(mysamples_array)


                    

def percentil_iterator(percentils_list, myfiles, number_genes, number_samples, uppercent, norm_accross_sample, outfile2, outfile3):
    
    outlier_matrix = np.zeros((number_genes, number_samples))
    
    for perc in percentils_list:
        # Process the expression matrix with a particular percentile
        thisGeneArray, list_of_samples, mysamples_array = create_mysamples_atpercentile(myfiles, perc, normal_samples)
        # Produce a list, and rank of outliers at that particular percentile, and add the rank to the outlier_matrix
        outlier_matrix = outliers_atpercentile(list_of_samples, uppercent,thisGeneArray.celllines, thisGeneArray.genelist, \
                                              outlier_matrix, norm_accross_sample)
    

    noutliergenes=10
    sample_outliers_dict=defaultdict(list)
    sample_outliers_dict_exp=defaultdict(list)
    sample_outliers_dict_med=defaultdict(list)
    sample_outliers_dict_std=defaultdict(list)

    
    indicators = np.array(range(number_genes))
    set_of_genes=set()
    genearray = np.array(thisGeneArray.genelist)
    folderpath='/home/oabalbin/projects/networks/Kinases_lists/output/'
    
    for j in range(number_samples):
        
        values = np.argsort(outlier_matrix[:,j])
        geneIndicators = indicators[values]
        geneIndicators = geneIndicators[-noutliergenes:]
        
        for gen in geneIndicators[::-1]:
            genename = thisGeneArray.genelist[gen]
            sample_outliers_dict[thisGeneArray.celllines[j]].append( [thisGeneArray.genelist[gen], outlier_matrix[gen,j]/len(percentils_list)])
            sample_outliers_dict_exp[thisGeneArray.celllines[j]].append( [thisGeneArray.genelist[gen], thisGeneArray.sampleCrudeExpValues[gen,j]] )
            sample_outliers_dict_med[thisGeneArray.celllines[j]].append( [thisGeneArray.genelist[gen], thisGeneArray.genemed[genename]] )
            sample_outliers_dict_std[thisGeneArray.celllines[j]].append( [thisGeneArray.genelist[gen], thisGeneArray.genesdt[genename]] )
            set_of_genes.add(thisGeneArray.genelist[gen])
            
            
        print thisGeneArray.celllines[j],'_sco',", ", ",".join(map(str,sum(sample_outliers_dict[thisGeneArray.celllines[j]],[]))).replace(",","\t")
        print thisGeneArray.celllines[j],'_raw',", ", ",".join(map(str,sum(sample_outliers_dict_exp[thisGeneArray.celllines[j]],[]))).replace(",","\t")
        print thisGeneArray.celllines[j],'_med',", ", ",".join(map(str,sum(sample_outliers_dict_med[thisGeneArray.celllines[j]],[]))).replace(",","\t")
        print thisGeneArray.celllines[j],'_std',", ", ",".join(map(str,sum(sample_outliers_dict_std[thisGeneArray.celllines[j]],[]))).replace(",","\t")
        
        outfile3.write(thisGeneArray.celllines[j]+'_sco'+"\t"+",".join(map(str,sum(sample_outliers_dict[thisGeneArray.celllines[j]],[]))).replace(",","\t")+'\n')
        outfile3.write(thisGeneArray.celllines[j]+'_raw'+"\t"+",".join(map(str,sum(sample_outliers_dict_exp[thisGeneArray.celllines[j]],[]))).replace(",","\t")+'\n')
        
        plot_sample_centric_view(outlier_matrix[geneIndicators[::-1],j], [], [], thisGeneArray.celllines[j], list(genearray[geneIndicators[::-1]]),\
                                  folderpath)    
            
      
    report_outlier_matrix = np.zeros((len(set_of_genes),number_samples))
    set_of_genes = list(set_of_genes)
    for cell, vect in sample_outliers_dict.iteritems():
        j = thisGeneArray.celllines.index(cell)
        for pair in vect:
            i = set_of_genes.index(pair[0])
            report_outlier_matrix[i,j] = pair[1]
    
    outfile2.write("geneName"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
    write_array2file(outfile2, report_outlier_matrix, set_of_genes)
    
        
    return sample_outliers_dict
    



def create_myArray(myfile, median_threshold, normal_samples):
    '''
    It creates a sample object and calculates a z-score for the expression level of each gene
    in the sample. Z-score (expvalues - genemedian)/genesd
    Returns a geneArray and also a samplelist with the expression for each sample
    The input matrix is already lo2 normalized    
    '''
    gene_dictionary, sample_dictionary, gene_interval_dict = defaultdict(), defaultdict(), defaultdict()
    celllines =  deque()
    normal_samples_list=[]
    
    first=True

    for line in myfile:
        if line is None:
            continue 
        if line.startswith('#'):
            continue
        # Handeling headers ....
       
        if line.startswith('gene'):
            celllines = line.strip('\n').split('\t')[2:]
            for j,sp in enumerate(celllines):
                sample_dictionary[sp] = j
                print sp,j
                if sp in normal_samples:
                    normal_samples_list.append(j)
                    
            k=0
            continue
        
        # Reading the file.
        fields = line.strip('\n').split('\t') 
        gene_dictionary[k] = fields[0]
        gene_interval_dict[k] = fields[1]
        
        exp=np.zeros(len(fields[2:]))
        expmasked = np.ma.make_mask_none((len(fields[2:]),))
        ## This function can be modified such as you verify all this after reading the whole
        ## array or better, making sure that all elemest in the array are already numbers.
        ## This could be a second function of verification
        
        exp = np.array(map(float,fields[2:]))
        
        '''
        for i, e in enumerate(fields[2:]):
            if e=="--" :
                expmasked[i] = False
                #exp[i] = float(0.1)
                exp[i] = float(1.0)
                #exp[i] = np.log2(float(0.1))
            else:
                #exp[i] = float(e) + 0.1
                exp[i] = float(e) + 1.0
                #exp.append(float(i))
        '''
        
        expvalues = np.array([exp]) + 1.0
        expvalues = np.ma.array(expvalues,mask=[expmasked])
        
        # shift the values by one. In order to avoid the zero coverage values which will be inf when log is applied.
        #expvalues = np.array([exp])
        #Creating an array of expression organized by sample
        # This could be better organized into a numpy array of defined dimensions.
        if first==True:
            geneRawValues = np.array(expvalues)
            first=False
        else:
            geneRawValues = np.concatenate((geneRawValues,expvalues),axis=0)
        
        k+=1
    
    # Pre-process the matrix data in order to filter out low expressing genes.
    # Calculate the median expression value for each gene.
    gene2keep = filter_low_expressing_genes(geneRawValues,median_threshold)
    geneRawValues = geneRawValues[gene2keep,:]
    gene_dictionary = rebuild_gene_dictionary(gene_dictionary,gene2keep)
    gene_interval_dict = rebuild_gene_dictionary(gene_interval_dict,gene2keep)
    
    # Create an GeneArray object
    thisGeneArray = myArray(geneRawValues, gene_dictionary, sample_dictionary, normal_samples_list, gene_interval_dict)

    return thisGeneArray

def filter_low_expressing_genes(geneRawValues,threshold):
    """
    It filters out genes with very low expression values. 
    genes which max expression value  < threshold for all samples are removed.          
    """
    print geneRawValues.shape[0]
    gene_indicators = np.arange(geneRawValues.shape[0])
    exp_cutoff = np.ma.median(geneRawValues, axis=1)
    #median_cutoff = sp.stats.scoreatpercentile(gene_medians,threshold) 
    #exp_cutoff = np.ma.max(geneRawValues, axis=1)
    gene2keep = gene_indicators[exp_cutoff > threshold]
    print len(gene2keep)
    return gene2keep

def rebuild_gene_dictionary(genedict,genes2keep):
    """
    """
    new_gene_dictionary=defaultdict()
    for i in range(len(genes2keep)):
        new_gene_dictionary[i] = genedict[genes2keep[i]]
    
    return new_gene_dictionary
    

def create_my_samples_list(thisGeneArray):
    
    med_percentile=25
    median_percentile = sp.stats.scoreatpercentile(genemed.values(),med_percentile)
    
    samplelist=[]
    mysamples_array= np.transpose(np.zeros((len(genelist) ,len(celllines))))
     
    for i,cell in enumerate(celllines):
        cellind = celllines.index(cell)
        geneDic, gene_norm_covDic, sample_norm_covDic= defaultdict(),defaultdict(),defaultdict()
        
        thissample_raw_values = sampleRawValues[:,cellind]
        #thissample_raw_values = samplegene_norm_cov[:,cellind]
        geneDic[cell] = thissample_raw_values #sampleRawValues[:,cellind]
        gene_norm_covDic[cell] = samplegene_norm_cov[:,cellind]
        
        #normalization of gene expression within a particular sample. 
        gene_median = np.median(thissample_raw_values)
        gene_mad = np.median(np.abs(thissample_raw_values - gene_median))   
        gene_std = np.std(thissample_raw_values)
        gene_insample_normalization = (thissample_raw_values - gene_median)/gene_mad
        #gene_insample_normalization = (thissample_raw_values - gene_median)/gene_mad
        
                
        sample_norm_covDic[cell] = gene_insample_normalization
        mysamples_array[i,:] = gene_insample_normalization
        
        thissample = Sample(cell, geneDic, genelist, gene_norm_covDic, sample_norm_covDic)
        samplelist.append(thissample)
    
    return thisGeneArray, samplelist, np.transpose(mysamples_array)


def specificity_sensitivity(normal_samples, outlier_samples, all_samples):
    """
    It calculates the specificity and sentivity of the outlier calling
    """
    # Specificity
    len_outlier_samples, len_normal_samples, len_all_samples =  float(len(outlier_samples)), float(len(normal_samples)), float(len(all_samples))
    false_positives = len(normal_samples.intersection(set(outlier_samples)))
    true_negatives = len_normal_samples - false_positives
    specificity = true_negatives / (true_negatives + false_positives)
    
    #Sensitivity
    true_positives = len_outlier_samples - false_positives
    false_negatives =  len_all_samples - len_normal_samples
    sensitivity = true_positives / (true_positives + false_negatives)
    #Prevalence
    prevalence = len_outlier_samples / (len_all_samples - len_normal_samples)
    
    return specificity, sensitivity, prevalence



############################
'''
thisfile = open('/home/oabalbin/projects/networks/Kinases_lists/output/mat_expfactor_allgenes_norm_raw_jun.txt')
#thisfile = open('/home/oabalbin/projects/networks/Kinases_lists/output/mat_expfactor_allgenesBreastCelllines_norm_raw_jun8.txt')
outfile_outliers= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersBreast_Expression_jun_10_100_m1.5_copascore_allgenes_mean_0.25_insample.txt','w')
outfile_outliers_names= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersBreast_names_jun_10_All_m1.5_copascore_allgenes_mean_0.25_insample.txt','w')
outfile_outliers_list = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersBreast_list_jun_10_100Outlier_m1.5_copascore_allgenes_mean_0.25_insample.txt','w')
outfile_outliers_matrix = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_Matrix_jun_10_200Outlier_m1.5_copascore_allgenes_mean_quality_25_insample.txt','w')
'''
'''
thisfile = open('/home/oabalbin/projects/networks/Kinases_lists/output/mat_expfactor_allgenesProstateTissue_norm_raw_jun8.txt')
##thisfile = open('/home/oabalbin/projects/networks/Kinases_lists/output/mat_expfactor_allgenesProstateTissue_norm_normals_jun8.txt')
outfile_outliers= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_Expression_jun_o8_200_m1.5_copascore_allgene_mean_quality.txt','w')
outfile_outliers_names= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_names_jun_o8_All_m1.5_copascore_allgene_mean_quality.txt','w')
outfile_outliers_list = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_list_jun_o8_200Outlier_m1.5_copascore_allgene_mean_quality.txt','w')
'''
'''
#thisfile = open('/home/oabalbin/projects/networks/EZH2_lists/output/mat_expfactor_allTUsProstateTissue_norm_raw_jun8.txt')
##thisfile = open('/home/oabalbin/projects/networks/Kinases_lists/output/mat_expfactor_allgenesProstateTissue_norm_normals_jun8.txt')
thisfile = open('/home/oabalbin/projects/networks/EZH2_lists/output/mat_expfactor_allTranscriptsProstateTissue_norm_raw_jun8.txt')
outfile_outliers= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_Expression_jun_10_200_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
outfile_outliers_names= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_names_jun_10_All_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
outfile_outliers_list = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_listB_jun_10_200Outlier_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
outfile_outliers_matrix = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_Matrix_jun_10_200Outlier_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
'''

thisfile = open('/home/oabalbin/projects/networks/EZH2_lists/output/allTranscripts_prostate_cohort.exprmat.txt')
outfile_outliers= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_Expression_jun_14_200_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
outfile_outliers_names= open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_names_jun_14_All_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
outfile_outliers_list = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_listB_jun_14_200Outlier_m1.5_copascore_AllTranscripts_mean_quality.txt','w')
outfile_outliers_matrix = open('/home/oabalbin/projects/networks/Kinases_lists/output/copa_outliersProstateTissue_Matrix_jun_14_200Outlier_m1.5_copascore_AllTranscripts_mean_quality.txt','w')





####################

normal_samples_list = list_of_names(open('/home/oabalbin/projects/networks/EZH2_lists/prostate_normal_samples.txt'))
#normal_samples_list = list_of_names(open('/home/oabalbin/projects/networks/Kinases_lists/breast_benigncelllines_samples.txt'))


gene_outliersample_list=[]
copa_score_genelist=[]
mean_gene_copa_score_list=[]
#percentils_list = [50]
percentils_list = [80,85,90,95,98] # [75,80,85,
gene_outliersample_dict, copa_rank_dict, copa_score_dict = defaultdict(list), defaultdict(list), defaultdict(list)
method='mean' # median, mean, sum
outliers_cutoff=200 #100 #200
nout2report=20
maximum_expression_cutoff=0.0 # 1.5 median RPKM
expression_cutoff=1.0
base_expression_cutoff= 50  # In percentil. It is used to report the expression value of a particular gen in a sample. 
                            # The median usually can be used, however because in many cases median=1, It is maybe desireable to use for e.g the75%. 
sdv_factor=0.25             #0.25
not_report_isoforms=True   # Report gene outlier isoforms

########### Program
#thisGeneArray = create_myArray(thisfile,maximum_expression_cutoff, normal_samples_list)
#pickle.dump(thisGeneArray,open('/home/oabalbin/projects/networks/Kinases_lists/output/kinasesBreast_array.tmp','w'))
#pickle.dump(thisGeneArray,open('/home/oabalbin/projects/networks/EZH2_lists/output/AllTranscriptsProstateTissuev17_5c_array.tmp','w'))
thisGeneArray = pickle.load(open('/home/oabalbin/projects/networks/EZH2_lists/output/AllTranscriptsProstateTissuev17_5c_array.tmp'))
#thisGeneArray.log2_norm()
#thisGeneArray.center_scale_norm()
thisGeneArray.copa_norm_across_samples(expression_cutoff)
thisGeneArray.copa_norm_in_sample()
#outliers_cutoff = thisGeneArray.expVal.shape[0]

#sys.exit(0)

tga = ga.gene_annotation("localhost", "oabalbin", "oscar", "gene_annotation")
############ Run the outlier analysis for different percentiles
for percentil in percentils_list:
    gene_outliersample_dict, copa_rank_dict, copa_score_dict = thisGeneArray.get_the_outliers(percentil, gene_outliersample_dict, copa_rank_dict, copa_score_dict)

############ Summarize the results from the different percentiles
indicator=np.array(range(len(percentils_list)))
gene_copascore_values = np.empty(len(copa_rank_dict))

for i in  range(len(copa_rank_dict)):
    if method == 'median':
        mean_copa_rank, mean_copa_score = np.median(copa_rank_dict[i]), np.median(copa_score_dict[i])
    elif method == 'mean':
        mean_copa_rank, mean_copa_score = np.mean(copa_rank_dict[i]), np.mean(copa_score_dict[i])
    #elif method == 'sum':
        
    genename, mean_copa_ind = thisGeneArray.gendict[i], indicator[copa_score_dict[i]>=mean_copa_score][0]
    genelocation = thisGeneArray.genInterval[i]
    outlier_samples = gene_outliersample_dict[i][mean_copa_ind]
    # quantifying how many normal samples the outlier hits
    normal_samples = set(thisGeneArray.normal_samples)
    all_samples = thisGeneArray.sampledict.values()
    # Calculation of specificity, sensitivy, and prevalence
    specificity, sensitivity, prevalence = specificity_sensitivity(normal_samples, outlier_samples, all_samples)
    quality_outliers = [specificity, sensitivity, prevalence]
    
    # Make a temporal array with only the copa score for each gene
    gene_copascore_values[i] = mean_copa_score
    # weigth the mean_copa_score for the number of samples in which that gene have a value greater than mean_copa_score
    mean_copa_score = mean_copa_score*len(outlier_samples)
    # Filter outliers according to the level of variance in the normals
    sdt_normals = thisGeneArray.sdv_normals(i)
    sdt_tumors = thisGeneArray.sdv_tumors_mets(i)
    if sdt_normals > sdt_tumors*sdv_factor:
        continue
    
    # Make a list of outliers  
    mean_gene_copa_score_list.append( gene_copa(i,genename,mean_copa_rank,mean_copa_score, mean_copa_ind, outlier_samples, \
                                                quality_outliers, genelocation) )
    

#sorted list of outliers
final_list_genes = sorted(mean_gene_copa_score_list, key=operator.attrgetter('mean_copa_score'), reverse=True)


############# Identify all outliers for each sample
outliers_in_samples=np.zeros((thisGeneArray.expVal.shape[0],thisGeneArray.expVal.shape[1]))
outliers_in_samples_copascore = np.zeros((thisGeneArray.expVal.shape[0],thisGeneArray.expVal.shape[1]))
list_unique_outliers = defaultdict()
for genename in thisGeneArray.gendict.values():
    list_unique_outliers[genename]=0

list_outliers_report = deque()
i=0
gene_translation_key, list_unique_outliers = thisGeneArray.gene_hugo_names(tga)


for gene in final_list_genes: 
#for i in range(len(final_list_genes)):
    if i > outliers_cutoff:
        break
    #gene = final_list_genes[i]    
    # Create a list of outliers to print and export to file
    #genename = gene_translation_key[gene.name]
    
    if gene.name.startswith('NM'):
        genename = tga.translate_refflat2hugo([gene.name])[0]
    else:
        genename = gene.name
    
    thisoutlier = [ i, gene.indicator, genename, gene.location, gene.mean_copa_score, np.median(thisGeneArray.expVal[gene.indicator]), 
                   np.max(thisGeneArray.expVal[gene.indicator]), np.min(thisGeneArray.expVal[gene.indicator]), gene.quality[0], gene.quality[1],
                   gene.quality[2] ]
    
    if not_report_isoforms:
        if list_unique_outliers[genename] == 0:
            list_outliers_report.append(thisoutlier)
            list_unique_outliers[genename] = 1
            i+=1
        else:
            continue
    else:
        list_outliers_report.append(thisoutlier)
        i+=1
        
    ### Most important thing  of this function    
    # Get the samples in which this gene is an outlier
    samples = gene_outliersample_dict[gene.indicator][gene.mean_copa_ind]
    outliers_in_samples[gene.indicator,samples] = 1
    # Matrix with copa score for the outliers in each sample
    outliers_in_samples_copascore[gene.indicator,samples] = gene_copascore_values[gene.indicator]

### Print outlier list and outlier matrix. Sample Indexes change the sample dictionary to be index-sample
print_samples=[]
for sp, j in thisGeneArray.sampledict.iteritems():
    print_samples.insert(j,sp)

outfile_outliers_matrix.write("geneName"+'\t'+",".join(print_samples).replace(',','\t')+'\n')

for thisoutlier in list_outliers_report:
    print ",".join(map(str,thisoutlier)).replace(',','\t')
    outfile_outliers_list.write(",".join(map(str,thisoutlier)).replace(',','\t')+'\n')
    
    eline = list(outliers_in_samples_copascore[thisoutlier[1],:])
    outfile_outliers_matrix.write(thisoutlier[2]+'\t'+",".join(map(str,eline)).replace(',','\t')+'\n') 
    
    
############# Print sample : outliers gene in that sample
gene_names_array = []
for i,gen in thisGeneArray.gendict.iteritems():
    gene_names_array.insert(i,gen)

#sys.exit(0)
gene_names_array = np.array(gene_names_array)
gene_indicators = np.array(range(thisGeneArray.expVal.shape[0]))
#median_gene_expression = np.ma.median(thisGeneArray.expVal,axis=1)
basal_gene_expression = thisGeneArray.get_gene_percentiles(base_expression_cutoff)

for sample,j in thisGeneArray.sampledict.iteritems():
    # gets the gene expression weighte / median expression for that gene in the population
    gene_expression_raw = np.array(thisGeneArray.expVal[outliers_in_samples[:,j]==1,j])/basal_gene_expression[outliers_in_samples[:,j]==1]
    # gets the raw copa score (not the weighted by num of samples). Rank using the copa score
    gene_expression = gene_copascore_values[outliers_in_samples[:,j]==1]
    #gene_expression = np.array(thisGeneArray.expVal_copainsample[outliers_in_samples[:,j]==1,j])
    # sort outliers genes according to expression
    sorted_outliers = np.ma.argsort(gene_expression)
    # Get gene names    
    oulier_genes_insample = gene_names_array[outliers_in_samples[:,j]==1]
    outliers_gene_expr = [[str(gene), str(exp)] for gene, exp in zip(oulier_genes_insample[sorted_outliers[::-1]], gene_expression_raw[sorted_outliers[::-1]])]
    outliers_gene_expr = sum(outliers_gene_expr,[])
    print j, sample, outliers_gene_expr[:nout2report]
    
    ### temporal
    outliers_gene_expr = [[str(gene), str(exp)] for gene, exp in zip(oulier_genes_insample[sorted_outliers[::-1]], gene_expression[sorted_outliers[::-1]])]
    outliers_gene_expr = sum(outliers_gene_expr,[])
    print j, sample, outliers_gene_expr[:nout2report]
    
    outfile_outliers_names.write(str(j)+'\t'+sample+'\t'+",".join(outliers_gene_expr[:nout2report]).replace(',','\t')+'\n')
    outfile_outliers.write(str(j)+'\t'+sample+'\t'+",".join(map(str,gene_expression[sorted_outliers[::-1]])).replace(',','\t')+'\n')
    
    
    

sys.exit(0)

folderpath='/home/oabalbin/projects/networks/Kinases_lists/output/'
myfiles = open_files('/home/oabalbin/projects/networks/Kinases_lists/output/','mat_expfactor_allgenes_norm_raw_jun.txt')#mat_expfactor_allgenes_norm_raw.txt')   

normal_samples = ['H16N2','HBL-100','HMEC1','HMEC2','MCF10A','MCF10F','MCF12A','hTERT-HME1']
                   
thisGeneArray, list_of_samples, mysamples_array = create_mysamples(myfiles, normal_samples)     

uppercent = 95
norm_accross_sample=False
#percentils_list = [50,70,80,90,95]
percentils_list = [50]
number_genes =  len(thisGeneArray.genelist)
number_samples = len(thisGeneArray.celllines)
outfile_cOutliers = open('/home/oabalbin/projects/networks/Kinases_lists/output/outlier_matrix_common50-95_mad_jun3.txt','w')
outfile_spOutliers=open('/home/oabalbin/projects/networks/Kinases_lists/output/outlier_matrix_common50-95_mad_top10Outliers_jun3.txt','w')
percentil_iterator(percentils_list, myfiles, number_genes, number_samples, uppercent, norm_accross_sample,outfile_cOutliers, outfile_spOutliers)


sys.exit(0)
# Note that this function modifies the sample objects in sample list
# the function outliers create a list of candidates in the sample object.
norm_accross_sample=False
outlier_matrix = outliers(list_of_samples, uppercent,thisGeneArray.celllines, thisGeneArray.genelist, norm_accross_sample)

# print an array of gene name and sample name. Presence or absence of that particular outlier in the sample
outfile = open('/home/oabalbin/projects/networks/Kinases_lists/output/outlier_matrix_raw_mod2_750_mad.txt','w')
outfile.write("geneName"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
write_array2file(outfile, outlier_matrix, thisGeneArray.genelist)
## Writing the expression matrices based on different normalizations

outfile1 = open('/home/oabalbin/projects/networks/Kinases_lists/output/expression_matrix_750_mad.txt','w')
outfile1.write("geneName"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
write_dict2file(outfile1, thisGeneArray.geneRawValues)

outfile2 = open('/home/oabalbin/projects/networks/Kinases_lists/output/median_score_matrix_750_mad.txt','w')
outfile2.write("geneName"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
write_dict2file(outfile2, thisGeneArray.geneExpInSample)

outfile3 = open('/home/oabalbin/projects/networks/Kinases_lists/output/zscore_score_matrix_750_mad.txt','w')
outfile3.write("geneName"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
write_dict2file(outfile3, thisGeneArray.geneZscoreInSample)

# Print the matrix with the sample normalization.
outfile4 = open('/home/oabalbin/projects/networks/Kinases_lists/output/zscore_sample_matrix_750_mad.txt','w')
print_sample_normalized_matrix(list_of_samples,thisGeneArray.celllines,thisGeneArray.genelist,outfile4)
#plot_outliers(list_of_samples,thisGeneArray,folderpath)

####### Matrices with only the expression value for the genes present in factors.

genes_in_factor= list_of_names_dict(open('/home/oabalbin/projects/networks/Kinases_lists/factor1.txt'))
'''
genes_pairs = sorted(genes_in_factor0.items(), key=operator.itemgetter(1), reverse=True)
genes_in_factor = defaultdict()
for gp in genes_pairs:
    print gp
    genes_in_factor[gp[0]] = gp[1]

print genes_in_factor
'''
outfile5 = open('/home/oabalbin/projects/networks/Kinases_lists/output/zscore_score_matrix_factor1c_750_mad.txt','w')
outfile5.write("geneName"+'\t'+"load"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
write_dict2file2(outfile5, thisGeneArray.geneZscoreInSample,genes_in_factor)

outfile6 = open('/home/oabalbin/projects/networks/Kinases_lists/output/median_score_matrix_factor1c_750_mad.txt','w')
outfile6.write("geneName"+'\t'+"load"+'\t'+",".join(thisGeneArray.celllines).replace(',','\t')+'\n')
write_dict2file2(outfile6, thisGeneArray.geneExpInSample,genes_in_factor)

outfile7 = open('/home/oabalbin/projects/networks/Kinases_lists/output/zscore_sample_matrix_factor1c_750_mad.txt','w')
print_sample_normalized_matrix(list_of_samples,thisGeneArray.celllines,thisGeneArray.genelist,outfile7,genes_in_factor)
















#gene="2_ERBB2"
#sample_group_erbb2 = sample_group_by_outlier_gene(samplelist, gene)
#print sample_group_erbb2
#recurrent_set_of_genes(samplelist)

    
        
        
        
        
        
        
        
        
        

